Freeze-out Conditions in Heavy Ion Collisions from QCD Thermodynamics 
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We present a determination of chemical freeze-out conditions in heavy ion collisions based on ratios 
of cumulants of net electric charge fluctuations. These ratios can reliably be calculated in lattice 
QCD for a wide range of chemical potential values by using a next-to-leading order Taylor series 
expansion around the limit of vanishing baryon, electric charge and strangeness chemical potentials. 
From a computation of up to fourth order cumulants and charge correlations we first determine 
the strangeness and electric charge chemical potentials that characterize freeze-out conditions in 
a heavy ion collision and confirm that in the temperature range 150 MeV < T < 170 MeV the 
hadron resonance gas model provides good approximations for these parameters that agree with 
QCD calculations on the (5-15)% level. We then show that a comparison of lattice QCD results for 
ratios of up to third order cumulants of electric charge fluctuations with experimental results allows 
to extract the freeze-out baryon chemical potential and the freeze-out temperature. 

PACS numbers: 11.15.Ha, 12.38.Gc, 12.38Mh, 24.60-k 



43- 



> 
O 
(N 

O 

(N 



X 



1 ) Introduction: A central goal of experiments at the 
Relativistic Heavy Ion Collider (RHIC) jl| is the ex- 
ploration of the phase diagram of Quantum Chromody- 
namics (QCD) at non-zero temperature (T) and baryon 
chemical potential (p-b)- In particular, a systematic 
Beam Energy Scan (BES) is being performed at RHIC in 
order to search for evidence for or against the existence 
of the QCD critical point, a second order phase tran- 
sition point, that has been postulated to exist at non- 
vanishing baryon chemical potential in the T-[j,b phase 
diagram of QCD 0, It would be the endpoint of a 
line of first order phase transitions which then would 
exist for larger [Ib- The measurement of fluctuations 
of conserved charges, e.g. net baryon number, electric 
charge and strangeness [H-[6j], plays a crucial role Q in 
this search for critical behavior and the exploration of 
the QCD phase diagram in general. 

Fluctuations of conserved charges generated in a heavy 
ion collision experiment may reflect thermal conditions 
at the time where the expanding medium, created in 
these collisions, cooled down and diluted sufficiently so 
that hadrons form again. It may be questioned whether 
the thermal medium at this time is in equilibrium and 
whether hadronization of all species takes place at the 
same time. However, statistical hadronization mod- 
els, based on thermal hadron distributions given by 
the Hadron Resonance Gas (HRG) model, describe the 
hadronization process quite successfully [8[- Moreover, 
HRG model calculations of net baryon number fluctua- 
tions Q describe well experimental data on net proton 
fluctuations This seems to suggest that at the time 
of chemical freeze-out the system can be described by 
thermodynamics characterized by a temperature Tf and 
a baryon chemical potential fi B . 



The measurement of conserved net charge fluctuations 
can provide a sensitive probe for critical behavior in hot 
and dense nuclear matter only when these fluctuations 
are generated at a point in the QCD phase diagram, char- 
acterized by (Tf, n B ), that is close to the QCD transi- 
tion line and eventually is also close to the elusive critical 
point. Lattice QCD calculations provide some informa- 
tion on the location of the QCD transition line in the 
T-/iB plane at small values of the baryon chemical po- 
tential 13, The position of the freeze-out points in 
this phase diagram are usually determined by compar- 
ing experimental data on multiplicities of various hadron 
species with the HRG model calculation 0, [l2[ . In or- 
der to put these parameters on a firm basis and com- 
pare them with the QCD transition line it is desirable 
to extract the freeze-out parameters by comparing ex- 
perimental data with a QCD calculation. This requires 
observables which arc experimentally accessible and can 
also reliably be calculated in QCD. The fluctuations of 
conserved charges and their higher order cumulants form 
such a set of observables. While net baryon number fluc- 
tuations are experimentally accessible only through mea- 
surements of net proton number fluctuations which 
may cause some difficulties [l3, 14 1, electric charge fluc- 
tuations may be easier to analyze. We thus will focus on 
the latter. As an intermediate step one should also verify 
to what extent HRG model calculations and QCD calcu- 
lations of freeze-out parameters yield consistent results 
when applied to the same set of thermal observables. 

We present here a calculation of ratios of cumulants 
of net electric charge and net baryon number fluctua- 
tions that can be formed from the first three cumulants. 
They are related to mean (Mx) , variance (o~ x ) and 
skewness (Sx) of the corresponding charge distributions, 
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X = B, Q, S for baryon number, electric charge and 
strangeness, respectively. These ratios can be compared 
to HRG model calculations and will be used to extract 
the freeze-out temperature and baryon chemical poten- 
tial from corresponding experimental measurements. 
2) Strangeness and electric charge chemical potentials: 
In order to get access to Tf and \J B we need to fix the 
electric charge (/xq) and strangeness (/«s) chemical po- 
tentials that characterize a thermal system created in a 
heavy ion collision. They are determined by assuming 
that the thermal sub-volume, probed by measuring fluc- 
tuations in a certain rapidity and transverse momentum 
window, reflects the net strangeness content and electric 
charge to baryon number ratios of the incident nuclei, 

M s = , M Q = rM B , (1) 

where M x = {VT^dln Z{n,T)/dfi x is the expec- 
tation value of the density of net charge X, \i = 
(pB, fJ-Q, fJ-s) summarizes the three charge chemical po- 
tentials and fix = (J-x/T. At any value of (T,/ib) the 
chemical potentials (/xq, fj,s) that fulfill these constraints 
can be evaluated in QCD. We perform a Taylor expansion 
of the densities Mx in terms of the three chemical poten- 
tials and calculate the expansion coefficients of this series 
using the lattice rcgularization scheme. This involves the 
numerical calculation of generalized susceptibilities 
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at a = 0. The calculation!]] of Y^? S becomes computa- 
tionally demanding at higher orders, but is done in lattice 
QCD with steadily increasing precision since many years 
[l5|. In particular, recent calculations performed with 
the Highly Improved Staggered Quark (HISQ) action flij 
and the stout action (l7j provide continuum extrapolated 
results for all diagonal (x*) and off-diagonal {Xii') sus- 
ceptibilities that are needed to determine (Aq, As) to 
leading order in As- 

Let us write the next-to-leading order (NLO) expan- 
sion of fiQ and /xs as 

Aq = Qi P-b + 93 As i As = si As + s 3 As ■ (3) 

Expanding the densities Mx up to third order in the 
chemical potentials we can fulfill the constraints speci- 
fied in Eq. [1] at NLO. This provides four equations to 
determine the four parameters (si, S3, q\, q^). In lead- 
ing order (LO) one obtains, 
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The NLO expressions are lengthy but can be derived eas- 
ily [2^ . We evaluated the leading order expressions in the 
temperature interval 150 MeV < T < 250 MeV for three 
different values of the lattice cut-off (a) corresponding to 
lattices with temporal extent N T = 1/aT = 6, 8 and 12. 
All calculations have been performed within an 0(a 2 ) im- 
proved gauge and staggered fermion (HISQ) discretiza- 
tion scheme [lj| for (2+l)-flavor QCD. The strange quark 
mass has been tuned to its physical value and the light 
to strange quark mass ratio is fixed to mi/m s = 1/20, 
which leads to a lightest Goldstone pion mass of about 
160 MeV. In the calculation of leading order results for 
As and Aq w e make use of data obtained by the HotQCD 
collaboration [l(|. On the 24 3 x 6, and 32 3 x 8 lattices we 
extended these calculations in the temperature interval 
150 MeV <T< 175 MeV to 30,000 molecular dynamics 
time units, saving configurations after every 10 , and by 
increasing the number of random vectors used to evaluate 
the susceptibilities to 1500 per gauge field configuration. 

We will in the following restrict our discussion to the 
case, r = 0.4, which approximates well the situation met 
in Au-Au as well as Pb-Pb collisions. The leading or- 
der expansion coefficients for Aq and As are shown in 
the top panels of Fig. [1] left and middle. Using spline 
interpolations of numerical results obtained for three dif- 
ferent lattice sizes, we performed extrapolations to the 
continuum limit using an ansatz linear in 1/N%. We have 
checked that no statistically significant differences occur 
by including an additional 1/N% correction. The result- 
ing extrapolations arc shown as bands in these panels. 

In order to check the importance of NLO corrections 
we have calculated S3 and 53 on lattices with temporal 
extent N T = 6 and 8. The results, expressed in units 
of the leading order terms, are also shown in Fig. [1] 
It is obvious from this figure that NLO corrections in- 
deed are small. They are negligible in the high tem- 
perature region and are below 10% in the temperature 
interval relevant for the analysis of freeze-out conditions, 
i.e., T ~ (160 ± 10) MeV. In fact, in this temperature 
range the leading order lattice QCD results deviate from 
HRG model calculations expanded to the same order by 
less than 15%. The next to leading order corrections 
start to become smaller than the HRG model values 
for T>160 MeV. This further reduces the importance 
of NLO corrections. In this respect we note that in the 
HRG model the NLO expansion reproduces the full HRG 
result for Aq and As to better than 1.0% for all values of 
Hb/T < 1.3. Altogether, we thus expect that the NLO- 
truncatcd QCD expansion is a good approximation to the 
complete QCD results for Aq and As for Ms < 200 MeV. 

In order to address further systematic errors we note 
that in a lattice QCD study as ours which utilizes the 
staggered discretization scheme, the biggest cut-off ef- 
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FIG. 1. The leading and next-to-leading order expansion coefficients of the strangeness (left) and the negative of the electric 
charge chemical potentials (middle) versus temperature for r — 0.4. For si and qi the LO-bands show results for the continuum 
extrapolation. For S3 and qz we give an estimate for continuum results (NLO bands) based on spline interpolations of the 
N T — 8 data. Dashed lines at low temperature are from the HRG model and at high temperature from a massless, 3-fiavor quark 
gas. The right hand panel shows NLO results for us/fiB and \iqI\ib as function of [is for three values of the temperature. 



fects at non-vanishing lattice spacing are generally due 
to so-called taste violations. These give rise to a dis- 
torted hadron spectrum but mainly influence the pion 
sector [lj| . Correspondingly the electric charge suscepti- 
bilities will be the ones most sensitive to discretization ef- 
fects while the baryon and strangeness sectors are largely 
unaffected. At leading order these discretization effects 
have been eliminated by taking the continuum limit. At 
NLO taste violation effects show up in the electric charge 
sector, Fig. 1 (middle). However, as the corrections them- 
selves are already small, we expect their influence to be 
small. Furthermore, the taste violations can be mod- 
elled within the HRG model by replacing the pion mass 
with the average, root mean-square pion mass [16[ of our 
lattice calculations. Results obtained with this modi- 
fied spectrum suggest that taste violation effects are in- 
deed negligible in the NLO calculation of the strangeness 
chemical potential and lead to at most 5% systematic 
errors in fiq for fig < 200 MeV. 

Additional systematic errors arise from the fact that we 
perform calculations with degenerate light quark masses 
m u = m^, as is usually done in current lattice QCD cal- 
culations. As a consequence not all susceptibilities xfjk S 
are independent; actually there arc two constraints in 
leading order ( X f = 2 X f« - vfi S , x f = 2j$f- X f 1 flr ) 
and six constraints in next to leading order [20j. These 
constraints of course, do not hold in the HRG model and 
may be considered as an additional source for the distor- 
tion of the hadron spectrum. Imposing these constraints 
by hand in the HRG model calculations we find that qi 
and can change by up to 3% while modifications of Si 
and s 3 are below the 1% level. This suggests that even 
after extrapolating to the continuum limit, the current 
lattice QCD calculations of hq/hb do have an inherent 
systematic error of about 3%. 



Our results for the strangeness and electric charge 
chemical potentials at NLO as function of fig and T are 
shown in Fig. QJright). While hs/^b varies between 0.2 
and 0.3 in the interval 150 MeV <T< 170 MeV, the ab- 
solute value of hq I /j,b is an order of magnitude smaller. 
Both ratios are almost constant for \lb < 200 MeV, which 
is consistent with HRG model calculations. 
3) Ratios of cumulants of net charge fluctuations: We 
now are prepared to evaluate cumulants of net charge 
fluctuations at as function of T and /jb at non- vanishing 
values for ^5 and /iq that obey the constraints appropri- 
ate for thermal conditions met in a heavy ion collision, 
i.e., Eq. [TJ Of particular interest are ratios of cumulants, 
JJ^ m = Xn,u,/Xm,u,> w hich to a large extent eliminate the 
dependence of cumulants on the freeze-out volume. Ra- 
tios with n + m even are non-zero for /.t = 0, while the 
odd-even ratios are in leading order proportional fis and 
thus vanish for /.t = 0. Ratios with n + m even or odd 
thus provide complementary information on Tf and pJ B . 
We will concentrate here on the simplest such ratios, 
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with X = B, Q. These ratios can be calculated in QCD 
as well as in the HRG model Q, and eventually can be 
compared to experimental data in order to determine Tf 
and n B . We evaluated them up to 0(p, B ) in a Taylor 
expansion for R* 2 and to leading order for . 

Let us first discuss the odd-even ratios Rf 2 - Using 
Xii X = X2 the LO coefficients can be written as 
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FIG. 2. The leading (top) and next-to-leading (bottom) order 
expansion coefficients of the ratio of first to second order cu- 
mulants of net electric charge fluctuations versus temperature 
for r — 0.4. The bands and lines are as in Fig. rjjleft). 

with X = B, Q. They have been evaluated on lattices 
with temporal extent N T = 6,8 and 12 and have been ex- 
trapolated to the continuum limit in the same way as for 
qi and si. The LO ratio Rf^ 1 , evaluated for three differ- 
ent values of the lattice cut-off (data points) and the re- 
sulting continuum extrapolation (LO band) are shown in 
Fig. EJtop). In Fig. [^bottom) we show the NLO correc- 
tions which have been evaluated on lattices with temporal 
extent N T = 6 and 8. The NLO corrections to the ratio 
of electric charge cumulants are below 10%, which makes 
the leading order result a good approximation for a large 
range of fiB- Systematic errors arising from the trunca- 
tion of the Taylor series for R,f 2 at next-to-leading order 
may again be estimated by comparing the full result in 
the HRG model calculation with the corresponding trun- 
cated results. Here we find for T = (160 ± 10) MeV and 
Hb/T < 1.3 that the difference is less than 1.0%. More- 
over, we estimated that taste violation effects in the NLO 
calculation lead to systematic errors that are at most 5% 
and thus will be negligible in Rf 2 . Taylor series truncated 
at NLO are thus expected to give a good approximation 
to the full result for a wide range of baryon chemical po- 
tentials. Similar results hold for the ratio Rf 2 , although 
NLO corrections are larger in this case. 

4-) Determination of freeze- out baryon chemical poten- 
tial and temperature: Obviously the ratio R,f 2 shows 
a strong sensitivity on [Ib but varies little with T for 
T ~ (160 ± 10) MeV. We show this ratio, evaluated in 
this temperature interval in a NLO Taylor expansion, 
in Fig. |3{left) as function of fis- For the determina- 
tion of (Tf,n B ) a second, complimentary information 
is needed. To this end we use the ratio R^ lt which is 
strongly dependent on T but receives corrections only at 
0{ji 2 B ). The leading order result for this ratio is shown 
in Fig. [3j middle). Apparently this ratio shows a charac- 
teristic temperature dependence for T>155 MeV that is 
quite different from that of HRG model calculations. The 



NLO correction to this ratio vanishes in the high temper- 
ature limit and at low T the HRG model also suggests 
small corrections. In fact, in the HRG model the LO 
contributions to E§y differ by less than 2% from the ex- 
act results on the freeze-out curve for fis < 200 MeV. In 
the transition region a preliminary 6 th order calculation 
at T = 162 MeV [2(| suggests that the /i B correction in 
units of the LO term is —0.03(10). Based on these esti- 
mates in the three T regions we expect the NLO correc- 
tions to be of the order of 10% for the whole temperature 
range. In Fig. [^middle) we show the spline interpolation 
for the N T = 8 data as a band and added on top of this 
a band that estimates the effect of a 10% contribution of 
the NLO correction. The ratio thus seems to be well 
suited for a determination of the freeze-out temperature. 

We now are in the position to extract fi B and Tf from 
Rf 2 and which eventually will be measured in the 
BES at RHIC. A large value for R^, i.e. R^ ~ 2 would 
suggest a low freeze-out temperature T<155 MeV, while 
a value i?^ ~ 1 would suggest a large freeze-out tem- 
perature, T ~ 170 MeV. A value of R^ ~ 1.5 would 
correspond to T ~ 160 MeV. 

A measurement of R^ thus suffices to determine the 
freeze-out temperature. In the HRG model parametriza- 
tion of the freeze-out curve [13] the favorite value for Tf 
in the beam energy range 200 GcV > v / s]vjv > 39 GeV 
indeed varies by less than 2 MeV and is about 165 MeV. 
At this temperature the values for i?^ calculated in the 
HRG model and in QCD differ quite a bit, as is obvious 
from Fig. El While R^ ~ 2 in the HRG model, one finds 
i$! ~ 1.2 in QCD at T = 165 MeV. Values close to the 
HRG value are compatible with QCD calculations only 
for T<157 MeV. We thus expect to cither find freeze-out 
temperatures that are about 5% below HRG model re- 
sults or values for i?^ that are significantly smaller than 
the HRG value. A measurement of this cumulant ratio 
at RHIC thus will allow to determine Tf and probe the 
consistency with HRG model predictions. 

For any of these temperature values a comparison of 
an experimental value for Rf 2 with Fig. (3Jlcft) will allow 
to determine \x B . To be specific let us discuss the results 
obtained at T = 160 MeV. Here we find 

R%{T = 160 MeV) = 0.102(5)As + 0.002(l)/t| . (8) 

For a value of the freeze-out temperature close to 
T = 160 MeV we thus expect to find n f B = (20-30) MeV, 
if Rf 2 lies in the range 0.012-0.020, fi f B = (50-70) MeV 
for 0.032 < R% < 0.045 and fi } B = (80 - 120) MeV for 
0.05 < i??2 < °- 08 - These parameter ranges are expected 
[H, [TH, [HI to cover the regions relevant for RHIC beam 
energies ^/sWn = 200 GcV, 62.4 GcV and 39 GeV, re- 
spectively. As is evident from Fig. [3^1cft) the values for 
H B will shift to smaller (larger) values when Tf turns 
out to be larger (smaller) than 160 MeV. A more refined 
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FIG. 3. The ratios Rf 2 versus ^lb/T (left) for three values of the temperature and versus temperature for /j,b = (middle). 
The wider band on the data set for N T = 8 (middle) shows an estimate of the magnitude of NLO corrections. The right hand 
panel shows the NLO result for the ratio of ratios of net electric charge and baryon number fluctuations, respectively. 



analysis of (X/,/Xg) will become possible, once the ratios 
R®2 and have been measured experimentally. 
5) Conclusions: We have shown that the first three cu- 
mulants of net electric charge fluctuations are well suited 
for a determination of freeze-out parameters in a heavy 
ion collision. Although the ratios Rf 2 and i?^ are suffi- 
cient to determine T/ and mu B , it will clearly be advan- 
tageous to have several ratios, including cumulants of net 
baryon number fluctuations, at hand that will allow to 
probe the consistency of an equilibrium thermodynamic 
description of cumulant ratios at the time of freeze-out. 
In particular the ratio of ratios Rf 2 /Ri2 = fY^./Y?.. is 
also well determined in (lattice) QCD calculations [161 ) . 
In Fig. [3jright) we show the NLO result for this ratio of 
ratios in the temperature range T = (160 ± 10) MeV. Its 
measurement will, on the one hand, allow to probe our 
basic assumptions on constraining the electric charge and 
strangeness chemical potentials and, on the other hand, 
constrain possible differences in cumulant ratios of net 
proton and net baryon number fluctuations. Once the 
ratios of lower order cumulants have been used to fix the 
freeze-out parameters, the calculation of higher order cu- 
mulants is parameter free and provides unique observ- 
ables for the discussion of possible signatures for critical 
behavior along the freeze-out line. 
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